home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / byt1186b.arc / RUNGKUT.LBR / ORBIT.FOR < prev    next >
Text File  |  1986-04-11  |  2KB  |  65 lines
  1. $NOFLOATCALLS
  2. $NODEBUG
  3. $STORAGE:2
  4. c****  User's calling program
  5. c****  NEQN=4 so dimension of work array = 4*17 = 68
  6.  
  7.        implicit double precision (a-h,o-z)
  8.        dimension y(4),work(68),icom(4)
  9.        external orbit
  10.        common alfasq
  11.  
  12.        open(2,file=' ',status='new')
  13.        icom(1)=0
  14.        write(*,*) 'imeth=,tola=,tolr='
  15.        read(*,*) imeth,tola,tolr
  16.        ecc=0.25d0
  17.        alfa=3.141592653589d0/4.d0
  18.        alfasq=alfa*alfa
  19.        neqn=4
  20.        hstart=0.01d0
  21.        y(1)=1.d0-ecc
  22.        y(2)=0.d0
  23.        y(3)=0.d0
  24.        y(4)=alfa*dsqrt((1.d0+ecc)/(1.d0-ecc))
  25.        x0=0.d0
  26.        xb=0.d0
  27.        icom(2)=0
  28.        icom(3)=0
  29.        do 20 j=1,24
  30.       xa=xb
  31.       xb=0.5d0*dble(j)+x0
  32.       write(2,100)xa,y(1),y(2)
  33.       call runkut(xa,y,xb,neqn,tola,tolr,hstart,work,
  34.      &               imeth,ierror,icom,orbit)
  35.       if(ierror.GT.1) then
  36.         write(2,100)xb,y(1),y(2)
  37.         write(2,*) ' ERROR-Problem too stiff or is discontinous'
  38.         close(2)
  39.         stop
  40.       end if
  41. 20     continue
  42.        if(icom(4).GT.0)write(2,*) ' Severe round-off error possible'
  43.        stop
  44. 100    format(F5.1,3F15.9)
  45.        end
  46. c********************************************************************
  47. c****  User supplied subroutine that contains the system of
  48. c****  differential equations to be solved.
  49. c****  Notice that in this routine it is necessary to have an array
  50. c****  yprime(neqn)
  51.  
  52.        subroutine orbit (x,y,yprime,neqn)
  53.        implicit double precision (a-h,o-z)
  54.        dimension y(neqn),yprime(neqn)
  55.        common alfasq
  56.  
  57.        r=y(1)*y(1)+y(2)*y(2)
  58.        r=r*dsqrt(r)/alfasq
  59.        yprime(1)=y(3)
  60.        yprime(2)=y(4)
  61.        yprime(3)=-y(1)/r
  62.        yprime(4)=-y(2)/r
  63.        return
  64.        end
  65.